home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / MCHB.FOR < prev    next >
Text File  |  1985-11-29  |  10KB  |  295 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE MCHB
  5. C
  6. C        PURPOSE
  7. C           FOR A GIVEN POSITIVE-DEFINITE M BY M MATRIX A WITH SYMMETRIC
  8. C           BAND STRUCTURE AND - IF NECESSARY - A GIVEN GENERAL M BY N
  9. C           MATRIX R, THE FOLLOWING CALCULATIONS (DEPENDENT ON THE
  10. C           VALUE OF THE DECISION PARAMETER IOP) ARE PERFORMED
  11. C           (1) MATRIX A IS FACTORIZED (IF IOP IS NOT NEGATIVE), THAT
  12. C               MEANS BAND MATRIX TU WITH UPPER CODIAGONALS ONLY IS
  13. C               GENERATED ON THE LOCATIONS OF A SUCH THAT
  14. C               TRANSPOSE(TU)*TU=A.
  15. C           (2) MATRIX R IS MULTIPLIED ON THE LEFT BY INVERSE(TU)
  16. C               AND/OR INVERSE(TRANSPOSE(TU)) AND THE RESULT IS STORED
  17. C               IN THE LOCATIONS OF R.
  18. C           THIS SUBROUTINE ESPECIALLY CAN BE USED TO SOLVE THE SYSTEM
  19. C           OF SIMULTANEOUS LINEAR EQUATIONS A*X=R WITH POSITIVE-
  20. C           DEFINITE COEFFICIENT MATRIX A OF SYMMETRIC BAND STRUCTURE.
  21. C
  22. C        USAGE
  23. C           CALL MCHB (R,A,M,N,MUD,IOP,EPS,IER)
  24. C
  25. C        DESCRIPTION OF PARAMETERS
  26. C           R      - INPUT IN CASES IOP=-3,-2,-1,1,2,3  M BY N RIGHT
  27. C                          HAND SIDE MATRIX,
  28. C                          IN CASE IOP=0  IRRELEVANT.
  29. C                    OUTPUT IN CASES IOP=1,-1  INVERSE(A)*R,
  30. C                           IN CASES IOP=2,-2  INVERSE(TU)*R,
  31. C                           IN CASES IOP=3,-3  INVERSE(TRANSPOSE(TU))*R,
  32. C                           IN CASE  IOP=0     UNCHANGED.
  33. C           A      - INPUT IN CASES IOP=0,1,2,3 M BY M POSITIVE-DEFINITE
  34. C                          COEFFICIENT MATRIX OF SYMMETRIC BAND STRUC-
  35. C                          TURE STORED IN COMPRESSED FORM (SEE REMARKS),
  36. C                          IN CASES IOP=-1,-2,-3  M BY M BAND MATRIX TU
  37. C                          WITH UPPER CODIAGONALS ONLY, STORED IN
  38. C                          COMPRESSED FORM (SEE REMARKS).
  39. C                    OUTPUT IN ALL CASES  BAND MATRIX TU WITH UPPER
  40. C                           CODIAGONALS ONLY, STORED IN COMPRESSED FORM
  41. C                           (THAT MEANS UNCHANGED IF IOP=-1,-2,-3).
  42. C           M      - INPUT VALUE SPECIFYING THE NUMBER OF ROWS AND
  43. C                    COLUMNS OF A AND THE NUMBER OF ROWS OF R.
  44. C           N      - INPUT VALUE SPECIFYING THE NUMBER OF COLUMNS OF R
  45. C                    (IRRELEVANT IN CASE IOP=0).
  46. C           MUD    - INPUT VALUE SPECIFYING THE NUMBER OF UPPER
  47. C                    CODIAGONALS OF A.
  48. C           IOP    - ONE OF THE VALUES -3,-2,-1,0,1,2,3 GIVEN AS INPUT
  49. C                    AND USED AS DECISION PARAMETER.
  50. C           EPS    - INPUT VALUE USED AS RELATIVE TOLERANCE FOR TEST ON
  51. C                    LOSS OF SIGNIFICANT DIGITS.
  52. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  53. C                     IER=0  - NO ERROR,
  54. C                     IER=-1 - NO RESULT BECAUSE OF WRONG INPUT
  55. C                              PARAMETERS M,MUD,IOP (SEE REMARKS),
  56. C                              OR BECAUSE OF A NONPOSITIVE RADICAND AT
  57. C                              SOME FACTORIZATION STEP,
  58. C                              OR BECAUSE OF A ZERO DIAGONAL ELEMENT
  59. C                              AT SOME DIVISION STEP.
  60. C                     IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  61. C                              CANCE INDICATED AT FACTORIZATION STEP K+1
  62. C                              WHERE RADICAND WAS NO LONGER GREATER
  63. C                              THAN EPS*A(K+1,K+1).
  64. C
  65. C        REMARKS
  66. C           UPPER PART OF SYMMETRIC BAND MATRIX A CONSISTING OF MAIN
  67. C           DIAGONAL AND MUD UPPER CODIAGONALS (RESP. BAND MATRIX TU
  68. C           CONSISTING OF MAIN DIAGONAL AND MUD UPPER CODIAGONALS)
  69. C           IS ASSUMED TO BE STORED IN COMPRESSED FORM, I.E. ROWWISE
  70. C           IN TOTALLY NEEDED M+MUD*(2M-MUD-1)/2 SUCCESSIVE STORAGE
  71. C           LOCATIONS. ON RETURN UPPER BAND FACTOR TU (ON THE LOCATIONS
  72. C           OF A) IS STORED IN THE SAME WAY.
  73. C           RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
  74. C           IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN RESULT MATRIX
  75. C           INVERSE(A)*R OR INVERSE(TU)*R OR INVERSE(TRANSPOSE(TU))*R
  76. C           IS STORED COLUMNWISE TOO ON THE LOCATIONS OF R.
  77. C           INPUT PARAMETERS M, MUD, IOP SHOULD SATISFY THE FOLLOWING
  78. C           RESTRICTIONS     MUD NOT LESS THAN ZERO,
  79. C                            1+MUD NOT GREATER THAN M,
  80. C                            ABS(IOP) NOT GREATER THAN 3.
  81. C           NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
  82. C           RESTRICTIONS ARE NOT SATISFIED.
  83. C           THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
  84. C           PARAMETERS ARE SATISFIED, IF RADICANDS AT ALL FACTORIZATION
  85. C           STEPS ARE POSITIVE AND/OR IF ALL DIAGONAL ELEMENTS OF
  86. C           UPPER BAND FACTOR TU ARE NONZERO.
  87. C
  88. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  89. C           NONE
  90. C
  91. C        METHOD
  92. C           FACTORIZATION IS DONE USING CHOLESKY-S SQUARE-ROOT METHOD,
  93. C           WHICH GENERATES THE UPPER BAND MATRIX TU SUCH THAT
  94. C           TRANSPOSE(TU)*TU=A. TU IS RETURNED AS RESULT ON THE
  95. C           LOCATIONS OF A. FURTHER, DEPENDENT ON THE ACTUAL VALUE OF
  96. C           IOP, DIVISION OF R BY TRANSPOSE(TU) AND/OR TU IS PERFORMED
  97. C           AND THE RESULT IS RETURNED ON THE LOCATIONS OF R.
  98. C           FOR REFERENCE, SEE H. RUTISHAUSER, ALGORITHMUS 1 - LINEARES
  99. C           GLEICHUNGSSYSTEM MIT SYMMETRISCHER POSITIV-DEFINITER
  100. C           BANDMATRIX NACH CHOLESKY - , COMPUTING (ARCHIVES FOR
  101. C           ELECTRONIC COMPUTING), VOL.1, ISS.1 (1966), PP.77-78.
  102. C
  103. C     ..................................................................
  104. C
  105.       SUBROUTINE MCHB(R,A,M,N,MUD,IOP,EPS,IER)
  106. C
  107. C
  108.       DIMENSION R(1),A(1)
  109.       DOUBLE PRECISION TOL,SUM,PIV
  110. C
  111. C        TEST ON WRONG INPUT PARAMETERS
  112.       IF(IABS(IOP)-3)1,1,43
  113.     1 IF(MUD)43,2,2
  114.     2 MC=MUD+1
  115.       IF(M-MC)43,3,3
  116.     3 MR=M-MUD
  117.       IER=0
  118. C
  119. C        MC IS THE MAXIMUM NUMBER OF ELEMENTS IN THE ROWS OF ARRAY A
  120. C        MR IS THE INDEX OF THE LAST ROW IN ARRAY A WITH MC ELEMENTS
  121. C
  122. C     ******************************************************************
  123. C
  124. C        START FACTORIZATION OF MATRIX A
  125.       IF(IOP)24,4,4
  126.     4 IEND=0
  127.       LLDST=MUD
  128.       DO 23 K=1,M
  129.       IST=IEND+1
  130.       IEND=IST+MUD
  131.       J=K-MR
  132.       IF(J)6,6,5
  133.     5 IEND=IEND-J
  134.     6 IF(J-1)8,8,7
  135.     7 LLDST=LLDST-1
  136.     8 LMAX=MUD
  137.       J=MC-K
  138.       IF(J)10,10,9
  139.     9 LMAX=LMAX-J
  140.    10 ID=0
  141.       TOL=A(IST)*EPS
  142. C
  143. C        START FACTORIZATION-LOOP OVER K-TH ROW
  144.       DO 23 I=IST,IEND
  145.       SUM=0.D0
  146.       IF(LMAX)14,14,11
  147. C
  148. C        PREPARE INNER LOOP
  149.    11 LL=IST
  150.       LLD=LLDST
  151. C
  152. C        START INNER LOOP
  153.       DO 13 L=1,LMAX
  154.       LL=LL-LLD
  155.       LLL=LL+ID
  156.       SUM=SUM+A(LL)*A(LLL)
  157.       IF(LLD-MUD)12,13,13
  158.    12 LLD=LLD+1
  159.    13 CONTINUE
  160. C        END OF INNER LOOP
  161. C
  162. C        TRANSFORM ELEMENT A(I)
  163.    14 SUM=DBLE(A(I))-SUM
  164.       IF(I-IST)15,15,20
  165. C
  166. C        A(I) IS DIAGONAL ELEMENT. ERROR TEST.
  167.    15 IF(SUM)43,43,16
  168. C
  169. C        TEST ON LOSS OF SIGNIFICANT DIGITS AND WARNING
  170.    16 IF(SUM-TOL)17,17,19
  171.    17 IF(IER)18,18,19
  172.    18 IER=K-1
  173. C
  174. C        COMPUTATION OF PIVOT ELEMENT
  175.    19 PIV=DSQRT(SUM)
  176.       A(I)=PIV
  177.       PIV=1.D0/PIV
  178.       GO TO 21
  179. C
  180. C        A(I) IS NOT DIAGONAL ELEMENT
  181.    20 A(I)=SUM*PIV
  182. C
  183. C        UPDATE ID AND LMAX
  184.    21 ID=ID+1
  185.       IF(ID-J)23,23,22
  186.    22 LMAX=LMAX-1
  187.    23 CONTINUE
  188. C
  189. C        END OF FACTORIZATION-LOOP OVER K-TH ROW
  190. C        END OF FACTORIZATION OF MATRIX A
  191. C
  192. C     ******************************************************************
  193. C
  194. C        PREPARE MATRIX DIVISIONS
  195.       IF(IOP)24,44,24
  196.    24 ID=N*M
  197.       IEND=IABS(IOP)-2
  198.       IF(IEND)25,35,25
  199. C
  200. C     ******************************************************************
  201. C
  202. C        START DIVISION BY TRANSPOSE OF MATRIX TU (TU IS STORED IN
  203. C        LOCATIONS OF A)
  204.    25 IST=1
  205.       LMAX=0
  206.       J=-MR
  207.       LLDST=MUD
  208.       DO 34 K=1,M
  209.       PIV=A(IST)
  210.       IF(PIV)26,43,26
  211.    26 PIV=1.D0/PIV
  212. C
  213. C        START BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  214.       DO 30 I=K,ID,M
  215.       SUM=0.D0
  216.       IF(LMAX)30,30,27
  217. C
  218. C        PREPARE INNER LOOP
  219.    27 LL=IST
  220.       LLL=I
  221.       LLD=LLDST
  222. C
  223. C        START INNER LOOP
  224.       DO 29 L=1,LMAX
  225.       LL=LL-LLD
  226.       LLL=LLL-1
  227.       SUM=SUM+A(LL)*R(LLL)
  228.       IF(LLD-MUD)28,29,29
  229.    28 LLD=LLD+1
  230.    29 CONTINUE
  231. C        END OF INNER LOOP
  232. C
  233. C        TRANSFORM ELEMENT R(I)
  234.    30 R(I)=PIV*(DBLE(R(I))-SUM)
  235. C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  236. C
  237. C        UPDATE PARAMETERS LMAX, IST AND LLDST
  238.       IF(MC-K)32,32,31
  239.    31 LMAX=K
  240.    32 IST=IST+MC
  241.       J=J+1
  242.       IF(J)34,34,33
  243.    33 IST=IST-J
  244.       LLDST=LLDST-1
  245.    34 CONTINUE
  246. C
  247. C        END OF DIVISION BY TRANSPOSE OF MATRIX TU
  248. C
  249. C     ******************************************************************
  250. C
  251. C        START DIVISION BY MATRIX TU (TU IS STORED ON LOCATIONS OF A)
  252.       IF(IEND)35,35,44
  253.    35 IST=M+(MUD*(M+M-MC))/2+1
  254.       LMAX=0
  255.       K=M
  256.    36 IEND=IST-1
  257.       IST=IEND-LMAX
  258.       PIV=A(IST)
  259.       IF(PIV)37,43,37
  260.    37 PIV=1.D0/PIV
  261.       L=IST+1
  262. C
  263. C        START BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  264.       DO 40 I=K,ID,M
  265.       SUM=0.D0
  266.       IF(LMAX)40,40,38
  267.    38 LLL=I
  268. C
  269. C        START INNER LOOP
  270.       DO 39 LL=L,IEND
  271.       LLL=LLL+1
  272.    39 SUM=SUM+A(LL)*R(LLL)
  273. C        END OF INNER LOOP
  274. C
  275. C        TRANSFORM ELEMENT R(I)
  276.    40 R(I)=PIV*(DBLE(R(I))-SUM)
  277. C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
  278. C
  279. C
  280. C        UPDATE PARAMETERS LMAX AND K
  281.       IF(K-MR)42,42,41
  282.    41 LMAX=LMAX+1
  283.    42 K=K-1
  284.       IF(K)44,44,36
  285. C
  286. C        END OF DIVISION BY MATRIX TU
  287. C
  288. C     ******************************************************************
  289. C
  290. C        ERROR EXIT IN CASE OF WRONG INPUT PARAMETERS OR PIVOT ELEMENT
  291. C        LESS THAN OR EQUAL TO ZERO
  292.    43 IER=-1
  293.    44 RETURN
  294.       END
  295.